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Abstract. Conventional methods for reducing the pulsed thermographic responses of delaminations tend to overestimate the 
size of the flaw. Since the heat diffuses in the plane parallel to the surface, the resulting temperature profile over the flaw 
is larger than the flaw. A variational method is presented for reducing the thermographic data to produce an estimated size 
for the flaw that is much closer to the true size of the flaw. The size is determined from the spatial thermal response of the 
exterior surface above the flaw and a constraint on the length of the contour surrounding the flaw. The technique is applied to 
experimental data acquired on a flat bottom hole composite specimen. 
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INTRODUCTION 

An effective method of rapid inspection of large carbon fiber reinforced polymer (CFRP) composite structures is 
flash thermography. The increased utilization of composites in commercial and military aircraft increases the need 
for viable inspection techniques for ensuring aircraft safety and reliability. Thermography is effective for detection of 
fatigue and impact damage in CFRP composites[l, 2, 3, 4, 5, 6], A variety of data reduction techniques have been 
shown to improve flaw detection[7, 8]. To assess the significance of the damage, it is often essential to be able to 
characterize the location and size of the damage. Previous efforts to estimate the depth of the damage have focused 
on the time evolution of the thermal response [7, 9, 10]. Since the heat diffuses in a direction parallel to the surface as 
the thermal disturbance from a flaw diffuses to the surface, the size of the thermal response is larger than the size of 
the flaw. Methods for obtaining the correct flaw size from the thermal response have included estimating the size from 
the time evolution of the shape of the thermal response[ll, 12], from the gradient of the thermal response[13] and a 
deconvolution of the point spread function of the thermal response[14, 15]. 

The "deblurring" of the thermal response to compensate for the diffusion in the plane parallel to the surface is 
similar to the deblurring problem in image processing. An effective technique for deblurring in image processing 
is a variational method, where a blurred estimation for the true image is compared to the measured image. The 
variational method for image processing is well developed[16, 17]. The direct application of image processing 
deblurring techniques to deblurring of thermographic images is somewhat limited since the desired flaw shapes are 
binary images (ie. either a flaw exists at a particular location or it does not), rather than gray or color images as is 
the case for image processing. This paper therefore examines of variational approach where the thermal response is a 
"blurred" image of a binary image. The method is formulated in a manner similar to the image segmentation technique 
developed by Chan and Vese[18], 

A simple way to fabricate and characterize a flaw in a specimen is a flat bottom hole. While it does not accurately 
represent a delamination, a flat bottom hole specimen was used to test the data reduction technique. The solution of 
the variational problem requires iterative solutions to a partial differential equation. The variational method is easy to 
formulate if the blurring function is a Gaussian, therefore a simple model for the thermal response of a flat bottom hole 
based on the convolution of a Gaussian kernal is presented and compared to a finite element simulation of the thermal 
response in the following section. 
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FIGURE 1. Comparison of the approximate convolution solution to finite element results for slots 0.2 cm, 0.4 cm and 2.0 cm wide 
and 0.05 cm below the surface. Compared profiles are for 0.5 seconds following the instantaneous surface heating. The approximate 
convolution solutions based on Eq.(2) are labeled as ACS and the finite element solutions are labeled as FEM. The nominal response 
(no flaw) has been subtracted from the thermal responses leaving only the responses due to the slot. 


THERMAL RESPONSE OF FLAT BOTTOM HOLE 


Approximate Analytical Solution for Flat Bottom Hole 


The thermal response to an instantaneous flux at one surface of a single layer is well known. The first two terms of 
the series solution of early time response based on image sources are 

t( ' ) = S| (1+2 ' 4) ' (1) 

where K is the thermal conductivity, K is the thermal diffusivity, / is the thickness of the layer and F is the instantaneous 
flux applied to the front surface. A possible approximate solution for a flat bottom hole response can be found by 
placing two sources, one below and one above the front surface, at a distance equal to twice the depth of the top of the 
flat bottom hole with the shape of the flat bottom hole. The thermal response at the front surface is then given by 
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where u(x' ,y') is the shape of the top surface of the flat bottom hole. If u(x J ,/) is a plane (that is «(Y ,/) = 1), then 
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and Eq.(2) reduces to Eq.(l) the solution for a layer with a thickness equal to the depth of the plane. 

To assess the accuracy of this approximation, the thermal response as calculated from Eq.(2) was compared to the 
results of a finite element calculation for a long slot in a layer. The material properties of the layer were a thermal 
conductivity of 0.97 W/m/°K, specific heat of 1200 J/kg/°K and a density of 1600 kg/m 3 . These are values commonly 
used for the surface normal thermal properties of a composite. A comparison of the finite element results and the 
thermal profiles obtained for 0.05 cm deep slot in a plate 0.2 cm thick from Eq.(2) is shown in Fig. 1. As can be seen 



FIGURE 2. Calculated the values of b as determined from finite element simulations of a 2 cm wide slot at three different depths 
in a 0.2 cm thick plate. 


from that figure, for large flaws near the surface where the size of the slot is significantly larger than the effective 
diffusion length (2 (k t) 1 / 2 )), the maximum amplitude both the finite element simulation and the convolution are the 
same. Near the edges of the slot, the convolution approximation produces more spreading of the thermal response than 
the finite element simulation, suggesting the 4-Kt in the Gaussian convolution kernel is too large. 


Thermal Response Based on Finite Element Simulation of Slot in Plate 


An important attribute of Eq.(2) is the in-plane diffusion of the response is not dependent on the depth of the flat 
bottom hole, but only on the time of the measured response. This implies that two holes at different depths, but having 
the same shape, will have the same shape thermal responses at the surface at any given time. While the shapes are the 
same, the amplitudes of the responses will be different and depend on depth. Since Eq.(2) seems to give a reasonable 
approximate solution, finite element simulations were performed for 2 cm wide slots with top surfaces 0.05 cm, 0. 1 
cm and 0.15 cm below the surface in a 0.2 cm thick plate. From these results an empirical approximation was sought. 

The empirical approximation assumes the early time thermal response has a form similar to Eq.(2) 
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where h is a constant. The constant h can be determined by comparison of the calculation of the thermal response using 
Eq. (4) and the results of the finite element simulations where the finite element simulated thermal response of the slot 
is larger than the computational accuracy of the simulation. For a slot, the shape of thermal response as determined 
from Eq.(4) is 
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where xl and x2 are the positions of the left side and right side of the slots respectively. Since xl and x2 are known, 
the thermal responses calculated from the finite element simulation can be fitted to determine b. The calculated values 
for b for each time step of the finite element simulation are shown in Fig. 2. The value of b, as determined from fitting 
the finite element simulations, vary about 2.2 instead of a value of 4, the effective value of b in Eq.(2). 





FIGURE 3. Comparison of the approximate convolution solution to finite element results for slots of 0.2 cm, 0.4 cm and 2.0 
cm wide and 0.05 cm below the surface. Compared profiles are for 0.5 seconds following the instantaneous surface heating. The 
approximate convolution solutions based on Eq.(2) are labeled as ACS and the finite element solutions are labeled as FEM.The 
nominal solution has been subtracted from the total solutions leaving only the response due to the slot. 


Using a value of 2.2 for b in Eq.(4), the thermal responses were calculated for the slotted plate configurations 
discussed in the previous section and the results are compared with the finite element solutions in Fig. 3. As can be 
seen from a comparison of that figure and Fig. 1, Eq. (4) gives a better approximation of the thermal response than 
Eq.(2) when b is equal to 2.2. In particular it gives a better approximation near the edges of the slot. For all but the 
smallest slot, it is difficult to distinguish the finite element simulation results from the convolution results. Comparisons 
of the finite element simulations for deeper slots, (0. 1 cm and 0. 15 cm below surface) give similar results. 

In the following section this approximation for the thermal response is applied within a variational method formu- 
lation to estimate the size of the top of flat bottom holes. 


APPLICATION OF VARIATIONAL PRINCIPLE TO FLAW SIZING 

The integral equation representation of the variational principle in image processing is 

E (p{x,y)) = JJ F{p(x,y),d(x,y))dx dy + y g(p(x,y)), (6) 

where E(p(x,y )) is the energy to be minimized, p(x,y ) is the desired image, d(x,y ) is the acquired image, F(p,d) is 
the relationship between the acquired image and the desired image, g(p) is a constraint and y is a constant. 

For sizing of a flat bottom hole from a single thermographic image that occurs at time T, if one assumes the signal 
due to the flat bottom hole is given by D(u(x,y),x,y, t) where u(x,y) is the shape of the flat bottom hole, then a possible 
energy to be minimized is 

E{u(x,y)) = JJ^\{D(u(x,y),x,y,T) -T(x,y,x)) 2 + y 8(<j> (x,y))\V u(x,y)\]dx dy, (7) 


where T{x,y, t) is the measured thermal response at time T and the constraint is proportional to the circumference of 
the flaw. To simplify estimating the shape of the flat bottom hole, rather than solving for u(x,y ) directly, a method is 
used that is similar to that proposed for image segmentation by Chan and Vese[18], that is the outer boundary of the 


shape is represented by the zero level set of a Lipschitz function <j> (x.y). The shape of the top of the flat bottom hole is 
given by H((j)(x,y)), where H is the heavy side function. 

To simplify the analysis, at each time step the nominal response of the plate is subtracted from the simulation of 
the slot in the plate leaving the response from the slot. This is normally considered the contrast in the image, with the 
nominal response set to zero. The response of the flat bottom hole at a given time T is 
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where B = bKt and A(l, t), the amplitude of the response, is dependent only on the depth of the top of the flat bottom 
hole(Z) and time of the response and shows explicitly that the response shape is independent of the depth of the top of 
the flat bottom hole. To first order, A(l, t) is given by 


A(l, x)=2e~ 


(9) 


a value that is found when reducing the data and could be used to estimate the depth of the top of the flat bottom hole. 
Representing R((j)(x,y)) as G*H(<j>(x,y)), then Eq.(7) becomes 

E (0(x,y))= [[ [(G*H(Q(x,y)) -T(x,y, t)) 2 +y 8(4>{x,y))\VH(<j>(x,y))\\dx dy. (10) 

JJn 

Following the approach taken by Chan and Vese[18], an approximation is substituted for the heavy side function and 
the (j>(x,y ) is determined from associated Euler-Lagrange equation 

8e(4>(x,y))(G*(G*H e (4>(x,y))-T(x,y,T)) +ydiv( |v^(^’y)| ), = ° ( n ) 

where we have chosen a slightly different approximation for the heavy side function of 
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and it’s derivative 
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e = 1 /2 seems to yield relatively good results for the spatial resolution of the thermographic images of interest. 
Assuming an initial guess for the boundary of the upper surface of the flat bottom hole is a single large circle, and the 
initial guess for (j) (x.y) is a signed distance function from that boundary, the solution to Eq.( 1 1) is found numerically 
using an iterative finite difference scheme in a manner similar to [18] and [ 19]. 


SYSTEM FOR MEASUREMENT OF THERMAL RESPONSE 

The thermal response of the composite specimen was performed with a commercial flash thermographic measurement 
system[4, 7]. The front surface of the composite is heated with flash lamp. The flash duration has been measured to be 
approximately 0.008 second. Since the earliest thermal responses of interest occur approximately one tenth of a second 
after an instantaneous heat pulse, this is a good estimate of an impulse excitation. The thermal response was measured 
with a focal plane array infrared imager detector with an array size of 640x512 that operates in the 3-5 micrometer 
wavelength band. The approximate spatial resolution of the acquired thermography images was 0.04 cm. The imager 
digital output frame rate was 60 hertz and was connected to a real time digital image processor to acquire the output 
images. All of the data of interest was collected within the 2.5 seconds after the flash heating. 


COMPOSITE FLAT BOTTOM HOLE SPECIMEN 

The composite specimen was approximately 10 cm x 10 cm with at thickness of 0.212 cm. Nine approximate flat 
bottom holes were drilled in from the back side of the specimen with depths of approximately one fourth, one half 




(a) (b) 

FIGURE 4. Computed tomography and thermal response of composite specimen, a. Representation of the x-ray computed 
tomography data acquired on composite flat bottom hole specimen. The hole numbers used to identify the holes in the paper 
are superimposed on the image, b. The measured thermal response from specimen acquired 0.25 seconds after the flash heating. 


TABLE 1. Size and Depth of Holes in Composite Flat Bottom Hole 
Specimen Based on the Computed Tomography Data 


Hole 

Depth from 

Percent thickness 


Number 

front surface 

above hole 

Diameter 

1 

0.154 ±0.003 cm 

73% 

0.32 ± 0.01 cm 

2 

0.156 ±0.003 cm 

74% 

0.64 ± 0.01 cm 

3 

0.148 ±0.005 cm 

70% 

1.28 ±0.01 cm 

4 

0.105 ± 0.008 cm 

50% 

0.32 ± 0.01 cm 

5 

0.109 ±0.004 cm 

51% 

0.64 ± 0.01 cm 

6 

0.097 ± 0.008 cm 

46% 

1.29 ±0.01 cm 

7 

0.051 ± 0.019 cm 

24% 

0.32 ± 0.01 cm 

8 

0.053 ± 0.007 cm 

25% 

0.64 ± 0.01 cm 

9 

0.027 ±0.014 cm 

13% 

1.28 ±0.01 cm 


and three fourths of the thickness, and diameters of approximately 1.27 cm, 0.63 cm and 0.32 cm. A more accurate 
measurements of the depth and diameters were determined from x-ray computed tomography data shown in Fig. 4 and 
are given in Table 1 . For comparison, an early time thermal response is shown in the same figure. 

As can be seen from the computed tomography data in Fig. 4, the top of the flat bottom holes are not quite flat. From 
an examination of the standard deviations of the depth, the holes with a diameter of approximately 0.64 cm have the 
flattest tops and the largest holes have the least flat tops. The deeper the hole was drilled the less flat the top of the 
hole. This lack of flatness is reflected in the thermal response shown as Fig. 4. In the thermal image it is clear that the 
slight increase in material at the center of the largest hole results in a slightly cooler location at the center of the hole 
response. Measurements at later times have this same cool spot at the center of all the largest holes, however, it is not 
as significant in the shallower the holes. 
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FIGURE 5. The estimation of the diameter of the thermal response of flat bottom holes with tops approximately 0.1 cm beneath 
the surface based on the reduction of images acquired at different times. The linear fits of the early time values used to estimate the 
diameters of the flat bottom holes are also shown in the the figure as solid lines and the nominal diameters of the holes are shown 
as dashed lines. 


ESTIMATION OF DIAMETER OF FLAT BOTTOM HOLE FROM THERMOGRAPHIC 

RESPONSE 

Two methods were used to estimate the diameters of the flat bottom hole specimens. Both of the methods have the 
advantage of not requiring an estimation of the depth of the flaw. The first method was based on a technique presented 
by Almond and Lau[l 1], A region centered on the flat bottom hole 2.44 cm by 2.44 cm (61x61 pixels) of the thermal 
response corresponding to a single flat bottom hole was selected from the data. At each time step where the thermal 
response of the flat bottom hole was significantly greater than the noise in the image, the nominal and maximum 
response was determined from the region following Gaussian smoothing. For each time step the data was thresholded 
with a cutoff at the mean of the maximum and nominal response to give a binary image representing the number of 
points with amplitudes greater the half the maximum. The sum of the number of pixels above the threshold times the 
area of a pixel is taken to be the area of the response for a given time step. Assuming the flaw shape is circular, the 
estimated diameter of the circle is determined from this area. The diameters of the thermal response as measured for 
the flat bottom holes with the depths of approximately 0.1 cm are shown in Fig. 5. The early time dependence of the 
the diameter of the response tends to be linear. The diameter of the response at time zero extrapolated from a linear 
fit of the early time response is the estimated diameter of the flat bottom hole. In Fig. 5, the linear fits used to estimate 
the diameter of the holes are also shown. A difference between the time dependence of the diameter of the thermal 
response shown here and the results obtained by Almond and Lau is the two smaller holes have positive slopes, while 
in the paper by Almond and Lau, all of the slopes were negative. It should be noted that Almond and Lau assumed 
a flaw size that was relatively large compared to the effective diffusion length(2(icf) 1 / 2 ) which is approximately 0.3 
cm for this specimen. The estimated diameters and the associated errors based on the confidence interval of the fit are 
given in Table .2. 

The same data was analyzed with the variational method. For each region of the thermal response around the center 
of the flat bottom hole, the nominal response at each time step of the composite was subtracted from the data in the 
region, producing a "contrast" image. Each "contrast" image is input into Eq.(l 1) as T (x,y, t). The variational method 
gives an estimation of the shape of the hole for each time step. A typical estimated shape obtained for the 0.32 cm 
diameter hole, 0.105 cm below the front surface (hole #4) based on the data obtained at 1.0 seconds is shown in 




TABLE 2. Size and Depth of Holes in Composite Flat Bottom Hole Specimen Based on the Computed Tomography 
Data 


Hole 

Number 

Depth from 
front surface 

Diameter 
of hole 

Diameter estimated from full 
width at half maximum* 

Diameter estimated from 
variational analysis 

1 

0.154 ±0.003 cm 

0.32 ± 0.01 cm 

0.82 ± 0.09 cm 

0.4 ± 0.2 

2 

0.156 ±0.003 cm 

0.64 ± 0.01 cm 

0.58 ± 0.05 cm 

0.64 ± 0.02 

3 

0.148 ±0.005 cm 

1.28 ± 0.01 cm 

1.03 ±0.02 cm 

1.23 ± 0.007 

4 

0.105 ±0.008 cm 

0.32 ± 0.01 cm 

0.47 ± 0.02 cm 

0.34 ± 0.06 cm 

5 

0.109 ± 0.004 cm 

0.64 ± 0.01 cm 

0.581 ±0.007 cm 

0.62 ± 0.01 cm 

6 

0.097 ± 0.008 cm 

1.29 ± 0.01 cm 

1.144 ± 0.007 cm 

1.22 ±0.06 cm 

7 

0.051 ±0.019 cm 

0.32 ± 0.01 cm 

0.308 ± 0.003 cm 

0.37 ± 0.007 cm 

8 

0.053 ± 0.007 cm 

0.64 ± 0.01 cm 

0.58 ± 0.05 cm 

0.56 ± 0.02 cm 

9 

0.027 ±0.014 cm 

1.28 ± 0.01 cm 

1.168 ± 0.002 cm 

1.07 ± 0.03 cm 


* Using method for estimation of size from time response discussed in Almond and Lau[l 1] 


Fig. 6. The thermal data used to calculate this estimate, a calculation of the thermal response from the shape of the 
hole and the difference between the calculated and measured responses are also shown in the figure. As can be seen 
from the figure, the shape of the thermal response calculated from the estimated shape for the hole agrees with the 
measured response to within the "noise" floor of the data. The "noise" floor is not random noise, but is a result of 
the inhomogeneities in thermal response of the composite. Currently this process has not been optimized and requires 
approximately 20 seconds per image to estimate the shape. 

For each time step the shape of the flat bottom hole was estimated using the variational method. The shape is a 
binary image with ones for the flat bottom holes and zeros where there was no hole. The sum of the number of pixels 
in the shape of flaw times the area of a pixel is the estimate of the area of the area of the flat bottom hole based on 
a single time step. The estimated diameter of the flat bottom hole based on a single time step assumes the shape is 
circular. For each time step the estimated diameters of the flat bottom holes with the depths of approximately 0. 1 cm 
are shown in Fig.7. As can seen from the figure, all of the time steps give approximately the same diameters. Relatively 
flat slopes of the time dependencies were also seen for all of the other holes, except hole #1 (0.154 cm below surface 
and 0.32 cm in diameter). This is expected if the convolution given in Eq.(4) is a good approximation for the measured 
thermal response. To obtain the flat time dependency, the diffusivity parallel to the surface ( K xx ) was assumed to be 
approximately 0.025 cm 2 /sec or about 5 times the measured surface normal thermal diffusivity (K±).The estimated 
diameters and the associated errors based on the standard deviation of the measurements at each time step are also 
given in Table .2. 

The techniques both produce a systematic error when estimating the diameter of the 1.28 cm holes. For both 
techniques the estimated diameter is smaller than the true diameter. This is perhaps a result of the top of the hole 
not being quite flat. It was noted in the discussion of Fig. 4 that the lack of flatness results in lower values at the center 
of the holes than at the edge for early times in the thermal response, which would not be expected for an ideal fat 
bottom hole response. This lack of flatness may have resulted in the systematic error of 1 to 2 mm. 

For both techniques the estimates for the diameters of the mid size holes were most accurate. For both techniques 
the estimated diameters were within 0.5 mm of the diameters measured from the computed tomography data. For the 
two deeper flaws, where the earliest time steps with contrasts large enough to allow for a measurement, the shape of 
the thermal response does not have the appearance of flat bottom hole response, but more closely resembles a two 
dimensional Gaussian as can be seen in Fig. 6, the variational technique estimates were slightly better those based on 
the full width at half maximum. 

The variational estimates for the diameters of the smallest size holes were more accurate than the diameters 
estimated from the full width at half maximum. The full width at half maximum results in estimated diameters that 
are much larger than the size of the hole for the two deepest 0.32 diameter holes. This is perhaps partly a result 
of uncertainty in the estimation of the values of the maximums. However, it is also a result of the spreading of the 
thermal response in the plane parallel to the surface. The variational approach gives an estimate of the diameter that is 
in excellent agreement with the computed tomography estimate, while the deepest hole is in agreement, but has very 
large error bars. 




FIGURE 6. Comparison of measured and calculated thermal response of composite specimen, a. Thermal response at 1.0 second 
of flat bottom hole with a diameter of 0.319 cm, 0.99 cm below the front surface, b. Estimation of the shape of the top of the flat 
bottom hole based on variational method, c. Thermal response calculated based on shape estimated of the top of the flat bottom 
hole. d. Difference between the calculated thermal response and measured response. 
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FIGURE 7. The variational method estimation of the diameter of the flat bottom holes with tops approximately 0.1 cm beneath 
the surface based on the reduction of images acquired at different times. A dashed line indicates the nominal diameter of the hole. 


CONCLUSIONS 

A method has been presented for determining the size of a flaw from the thermographic response using a variational 
approach to reduce the data. The technique was applied to thermographic data acquired on a flat bottom hole composite 
specimen, with flaws at different depths and diameters. From the estimated shapes of the holes, the hole diameters 
were found to be in good agreement with measurements of the hole diameters from computed tomography data on 
the specimen. The results were also compared with the estimated diameters from the same thermographic data using a 




method described by Almond and Lau[ll]. There was good agreement between the estimated diameters from the two 
techniques for all but the smallest holes. 

One of the advantages of the variational technique is that it is easily applied to shapes other than circular flat 
bottom holes. The shape obtained from the technique should be independent of the time step of the thermal data 
reduced. It is possible to reformulate the technique to reduce multiple time steps at the same time, to improve the 
shape estimation and estimate the depth of the flaw. Future efforts will focus on development of this formulation and 
testing the application of the technique to sizing of delaminations, which are of more interest in practical applications. 
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